function [z, df, d2f] = FUNC5(a,b,c,d,e,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This version is checked by Jun Yu on Oct 31, 2023
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1

% checked
fun = @(x)exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
z = integral(fun,a,100);

% checked
da=-exp(-(a-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));

% checked
funb = @(x)1./(b).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
db=integral(funb,a,100);

% checked
func = @(x)1./(c.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
dc=integral(func,a,100);

% checked
fund = @(x)1./(d).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
dd=-integral(fund,a,100);

% checked
fune = @(x)(-exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.^2.*sqrt(2.*pi))+exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(sqrt(2.*pi)).*(-((x+0.5.*e.^2).*e.^2-(x+0.5.*e.^2).^2)./e.^4)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
de = integral(fune,a,100);

% checked
fung = @(x)exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(-exp(-(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2./(2.*f.^2))./(f.^2.*sqrt(2.*pi))+exp(-(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2./(2*f.^2))./(sqrt(2.*pi)).*(-((log(b.*c/(d.*(1+c)))-x+0.5.*f.^2).*f.^2-(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2)./f.^4));
dg = integral(fung,a,100);

df=[da db dc dd de dg];


% second order derivatives
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1
    % define the probablity density function
    function [f0] = desity(x,sigma_under)
      f0 = 1/(sigma_under*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2));
    end
    
    % call the probability density function
    f0_omega0bar = desity(a,e);
    f1_omega0bar = desity(log(b*c/(d*(1+c)))-a,f);
    
    % first order derivative of probability density function with respect
    % to sigma
    function [f0prime] = desityprime(x,sigma_under)
      f0prime = -1/(sigma_under^2*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))+ 1/(sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2))*(-(((x+0.5*sigma_under^2)*sigma_under^2 -(x+0.5*sigma_under^2)^2)/(sigma_under^4)));
    end
    
    % Second order derivative of probability density function
    function [f0prime2] = desityprime2(x,sigma_under)
      f0prime2 = 2/(sigma_under.^3.*sqrt(2.*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)) ...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*((((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)/(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-((2.*x.*sigma_under.^2+2.*sigma_under.^4-6.*sigma_under.^2.*(x+0.5.*sigma_under.^2)+4.*(x+0.5.*sigma_under.^2).^2)./(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-(((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)./(sigma_under.^4))).*(-(((x+0.5*sigma_under.^2).*sigma_under.^2 -(x+0.5*sigma_under.^2).^2)/(sigma_under.^3)));
    end
    f0prime_omega0bar = exp(-(a-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))*(-(a-(-0.5.*e.^2))./(e.^2));
    f1prime_omega0bar = exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2))./(f.^2)).*(-1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For da functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % get daa % checked
    daa = -f1prime_omega0bar*f0_omega0bar-f1_omega0bar*f0prime_omega0bar;
    
    % get dab % checked
    dab = -exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b*c/(d*(1+c)))-a-(-0.5.*f.^2))./(f.^2)).*(1./b).*f0_omega0bar;
    
    % get dac % checked
    dac = -exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2))./(f.^2)).*(1./(c.*(1+c))).*f0_omega0bar;
    
    % get dad % checked
    dad = -exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2))./(f.^2)).*(-1/d).*f0_omega0bar;
    
    % get dae % checked
    dae = -f1_omega0bar.*desityprime(a,e);
    % get dag % checked
    dag = -desityprime(log(b.*c./(d.*(1+c)))-a,f).*f0_omega0bar;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For db functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dba 
    dba = dab; % checked
    
    % get dbb % checked
    %Phsi = (-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2));
    funbb = @(x)(((-1./(b.^2)).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)) + (1./b).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./b).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./b).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(b.*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dbb = integral(funbb,a,100);
    
    % get dbc % checked
    funbc = @(x)(((1./b).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c))).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./b).*exp(-(log(b*c/(d*(1+c)))-a-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(c.*(1+c).*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dbc = integral(funbc,a,100);
    
    % get dbd % checked
    funbd = @(x)(((1./b).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./b).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(1./(d.*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dbd = integral(funbd,a,100);    
    
    % get dbe % checked, a problem found and corrected on Oct 31, 2023
    funbe = @(x)((1./b).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2))))...
                .*(-1./(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1/(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dbe = integral(funbe,a,100); 
    
    % get dbg % checked 
    funbg = @(x)(((1./b).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1/(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))))).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./b).*exp(-(log(b.*c/(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(2.*(log(b.*c./(d.*(1.+c)))-x)./(f.^3))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dbg = integral(funbg,a,100);  
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dc functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    % get dca 
    dca = dac; % checked
    
    % get dcb
    dcb = dbc; % checked
    
    % get dcc % checked
    funcc = @(x)(((-(1+2.*c)/((c.*(1+c)).^2)).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)) + (1./(c.*(1+c))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c))).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./(c.*(1+c))).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(c.*(1+c).*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dcc = integral(funcc,a,100);
    
    % get dcd % checked
    funcd = @(x)(((1./(c.*(1+c))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./(c.*(1+c))).*exp(-(log(b.*c/(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(1./(d.*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dcd = integral(funcd,a,100);     
    
    % get dce  % checked
    funce = @(x)((1./(c.*(1+c))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2))))...
                .*(-1./(e.^2.*sqrt(2.*pi))*exp(-(x+0.5.*e.^2).^2/(2.*e.^2))+ 1/(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dce = integral(funce,a,100); 
    
    % get dcg  % checked
    funcg = @(x)(((1./(c.*(1+c))).*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)/(f.^4))))).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              +(1./(c.*(1+c))).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(2.*(log(b.*c./(d.*(1.+c)))-x)./(f.^3))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    dcg = integral(funcg,a,100);      
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dd functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
    % get dda
    dda = dad;  % checked
    
    % get ddb
    ddb = dbd;  % checked
    
    % get ddc;
    ddc = dcd;  % checked
    
    % get ddd  % checked
    fundd = @(x)(((1./(d.^2)).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)) - (1./d).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              -(1./d).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(1/(d.*f.^2))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    ddd = integral(fundd,a,100);
    
    % get dde  % checked
    funde = @(x)((1./(d)).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2))))...
                .*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1/(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)/(e.^4))));
    dde = -integral(funde,a,100); 
    
    % get ddg  % checked
    fundg = @(x)(((-1./d).*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2)*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)/(f.^4))))).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2))...
              -(1./d).*exp(-(log(b.*c./(d.*(1+c)))-a-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(2.*(log(b.*c./(d.*(1.+c)))-x)./(f.^3))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)));
    ddg = integral(fundg,a,100);    
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For de functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dea
    dea = dae;  % checked
    
    % get deb
    deb = dbe;  % checked
    
    % get dec
    dec = dce;  % checked
    
    % get ded
    ded = dde;  % checked
    
    % get dee;  % checked
    funee = @(x)exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*desityprime2(x,e);
    dee = integral(funee,a,100); 

    % get deg  % checked
    funeg = @(x)(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4)))).*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    deg = integral(funeg,a,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dg functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dga
    dga = dag;  % checked
    
    % get dgb
    dgb = dbg;  % checked
    
    % get dgc
    dgc = dcg;  % checked
    
    % get dgd
    dgd = ddg;  % checked
    
    % get dge
    dge = deg;  % checked
    
    % get dgg  % checked
    % may be problematic
    fungg = @(x)exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi)).*desityprime2(log(b.*c./(d.*(1.+c)))-x,f);
    dgg = integral(fungg,a,100);   
    
d2f = [daa dab dac dad dae dag; 
       dba dbb dbc dbd dbe dbg;
       dca dcb dcc dcd dce dcg;
       dda ddb ddc ddd dde ddg;
       dea deb dec ded dee deg;
       dga dgb dgc dgd dge dgg];

end



